##########################################################################################

library(data.table)
library(optparse)
library(dplyr)
library(ggplot2)
library(ggsci)
library(ggpubr)

##########################################################################################

option_list <- list(
    make_option(c("--sample_list_file"), type = "character"),
    make_option(c("--rsem_file"), type = "character"),
    make_option(c("--gtf_file"), type = "character"),
    make_option(c("--gene"), type = "character"),
    make_option(c("--out_path"), type = "character")
)

if(1!=1){
    
    sample_list_file <- "~/20220915_gastric_multiple/rna_combine/analysis/config/tumor_normal.list"
    rsem_file <- "~/20220915_gastric_multiple/rna_combine/analysis/images/DiffGene/CombineCounts.FilterLowExpression-MergeMutiSample.varianceStabilizingTransformation.tsv"
    out_path <- "~/20220915_gastric_multiple/rna_combine/analysis/images/showGene"
    gtf_file <- "~/ref/GTF/gencode.v19.ensg_genename.txt"
    gene <- "FAP"

}

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

sample_list_file <- opt$sample_list_file
out_path <- opt$out_path
rsem_file <- opt$rsem_file
gtf_file <- opt$gtf_file
gene <- opt$gene

##########################################################################################

info <- data.frame(fread(sample_list_file))
dat_expression <- data.frame(fread(rsem_file))
dat_gtf <- data.frame(fread(gtf_file , header = F))

colnames(dat_gtf) <- c("gene_id" , "Hugo_Symbol")

##########################################################################################

igc_class <- c("IM + IGC + DGC" , "IM + IGC")
dgc_class <- c("IM + IGC + DGC" , "IM + DGC")

im_igc_sample <- unique(info[info$Type %in% igc_class , "ID"])
im_dgc_sample <- unique(info[info$Type %in% dgc_class , "ID"])

##########################################################################################
dat_expression <- merge( dat_expression , dat_gtf , by = "gene_id" )

##########################################################################################
plotTpm <- function(tmp_dat = tmp_dat , out_name = out_name , title = title){

    my_comparisons_1 <- list( 
        c(1, 2) , c(1, 3) , c(1, 4) , 
        c(2, 3) , c(2, 4) ,
        c(3, 4)  
        )

    plot <- ggplot( tmp_dat , aes( x = Class , y = Tpm , color = Class ) ) +
        geom_line( aes( group = Sample ) , size = 0.4 , color = "gray" ) +  ## 配对样本加线
        geom_boxplot(alpha =1 , outlier.size=0 , size = 0.9 , width = 0.6) +
        geom_jitter(position = position_jitter(0.17) , size = 1 , alpha = 0.7) +
        scale_fill_npg()+
        scale_color_npg()+
        xlab(NULL) +
        ylab("TMM")+
        theme_bw() +
        ggtitle(title) +
        stat_compare_means(comparisons = my_comparisons_1) +
        theme(panel.background = element_blank(),#设置背影为白色#清除网格线
            legend.position ='none',
            legend.title = element_blank() ,
            panel.grid.major=element_line(colour=NA),
            plot.title = element_text(size = 8,color="black",face='bold'),
            legend.text = element_text(size = 10,color="black",face='bold'),
            axis.text.y = element_text(size = 10,color="black",face='bold'),
            axis.title.x = element_text(size = 10,color="black",face='bold'),
            axis.title.y = element_text(size = 10,color="black",face='bold'),
            axis.ticks.x = element_blank(),
            axis.text.x = element_text(size = 10,color="black",face='bold') ,
            axis.line = element_line(size = 0.5)) 
    ggsave(file=out_name,plot=plot,width=4,height=5)

}

describeGene <- function(dat_expression = dat_expression , image_path = image_path ){

    tmp_exp <- subset( dat_expression , Hugo_Symbol == gene)

    sample <- sapply( strsplit(colnames(tmp_exp)[2:(ncol(tmp_exp)-1)] , "_") , "[" , 1)
    class <- sapply( strsplit(colnames(tmp_exp)[2:(ncol(tmp_exp)-1)] , "_") , "[" , 2)
    tpm <- as.numeric(tmp_exp[2:(ncol(tmp_exp)-1)])
    tmp_dat <- data.frame( Sample = sample , Class = class , Tpm = tpm )
    tmp_dat$Class <- factor( tmp_dat$Class , levels = c( "Normal" , "IM" , "IGC" , "DGC" ) , order = T )

    out_name <- paste0( image_path , "/" , gene , ".normalize.oneImage.pdf" )

    title <- paste0(
        "Gene : " , gene , "\n" 
    )
    ## 差异表达的P值
    plotTpm(tmp_dat = tmp_dat , out_name = out_name , title = title)

    out_name <- paste0( image_path , "/" , gene , ".normalize.tsv" )
    write.table(tmp_exp , out_name , row.names = F , sep = "\t" , quote = F)
}


##########################################################################################
## 分单个基因，看在一个人多个样本的改变情况
Hugo_Symbol <- gene
image_path <- out_path
dir.create(image_path , recursive = T)

describeGene(dat_expression = dat_expression , image_path = image_path )



